Limitations of Quantum Simulation Examined by Simulating a Pairing Hamiltonian 

using Nuclear Magnetic Resonance 
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Quantum simulation uses a well-known quantum system to predict the behavior of another quan- 
tum system. Certain limitations in this technique arise, however, when applied to specific problems, 
as we demonstrate with a theoretical and experimental study of an algorithm to find the low-lying 
spectrum of a Hamiltonian. While the number of elementary quantum gates required does scale 
polynomially with the size of the system, it increases inversely to the desired error bound e. Making 
such simulations robust to decoherence using fault-tolerance constructs requires an additional factor 
of ~ 1/e gates. These constraints are illustrated by using a three qubit nuclear magnetic resonance 
system to simulate a pairing Hamiltonian, following the algorithm proposed by Wu, Byrd, and Lidar 
[!]■ 



The unknown properties and dynamics of a given quan- 
tum system can often be studied by using a well-known 
and controllable quantum system to mimic the behav- 
ior of the original system. This technique of quantum 
simulation is one of the fundamental motivations for the 
study of quantum computation [2-4] , and is particularly 
of interest because a quantum simulation may be per- 
formed using space and time resources comparable to the 
original system. Such "efficient" scaling is dramatically 
better than the exponentially large resource requirements 
to simulate any general quantum system with a classical 
computer, as Feynman originally observed [2]. 

Recent work has continued to arouse great interest in 
quantum simulation, because it offers the possibility of 
solving computationally hard problems without requiring 
the resources necessary for algorithms such as factoring 
[5] and searching [6]. Experimental results have demon- 
strated simulations of a truncated oscillator and of a 
three-body interaction Hamiltonian, using a nuclear mag- 
netic resonance (NMR) quantum computer [7, 8], and 
explored various solid-state models on two qubit systems 
[9-12]. Interest has also extended to simulating complex 
condensed matter systems with quantum optical systems 
[13], demonstrated vividly by the observation of a su- 
pcrfiuid to Mott insulator transition in a Bose-Einstein 
condensate [14]. 

Often overlooked in the discussion of quantum simu- 
lations, however, is the question of desired precision (or 
error e) in the final measurement results. Current quan- 
tum simulation techniques generally scale poorly with de- 
sired precision; they demand an amount of space or time 
which increases as 1/e, broadly translating into a num- 
ber of quantum gates which grows exponentially with the 
desired number of bits in the final answer. Why is this 
scaling behavior so poor, and what is its physical origin? 

Consider as a specific example the problem of calculat- 
ing the energy gap A between the ground state \G) and 
the first excited state £a), of a Hamiltonian H. A can 
be found using the following steps: 1) map the Hilbert 



space of the system to be simulated to n qubits, 2) pre- 
pare the computer in the state [*/) = cg\G) +ce\E\), 
3) evolve under the Hamiltonian for times ti, 4) extract 
the phase difference as a function of time between the 
evolution of the ground and first excited state. 

Two methods for calculating the phase difference, and 
thus A, are as follows. The first method uses the phase 
estimation algorithm [15, 16]. This method relies on the 
quantum Fourier transform (QFT) and requires simulat- 
ing the Hamiltonian for times tk = 2 fc io, for integer k 
from to q. Since the input state is a superposition of 
\G) and \E\), the measured phase will either be Ecto or 
(Eq + A)to, where Eg is the ground state energy. One 
needs to run the algorithm on average 1 /c E times to get 
both values and thus measure A. 

The second method does not use the QFT, and in- 
stead simulates the Hamiltonian for times tk — kto, for 
integer k from to Q, and then measures any operator 
M such that (G\ M \Ei) ^ 0. Typically, any operator 
that does not commute with the Hamiltonian suffices. 
After calculating (M(tk)), one classically Fourier trans- 
forms (FT) over the averaged values yielding a spectrum 
(M(u))) with peaks at ±A and 0. 

For fixed precision, obtaining A up to error e for fixed 
e, both methods can be "efficient," in that the number of 
elementary steps (or quantum gates) required increases 
only polynomially with the number of qubits n, if the 
Hamiltonian can be efficiently simulated and the initial 
states efficiently prepared. A d-qubit Hamiltonian can be 
simulated with a number of gates of order 0(n d ) assum- 
ing two qubit interactions between any qubits [4]. If one 
assumes only nearest neighbor two qubit gates, it scales 
as 0(n d+1 ). Most physical systems of interest are de- 
scribed by two-body interactions which can be described 
by four qubit Hamiltonians. 

Consequently, the challenge of designing efficient quan- 
tum simulations is choosing a property that can be ef- 
ficiently extracted. However, no general measurement 
method is known which allows A to be measured effi- 
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ciently with respect to the precision using such quantum 
simulations. For error e, the number of digits of precision 
in the result is log(l/e), and both of the above methods 
require ~ 1/e steps (or gates) to obtain this precision. 
In contrast, an efficient algorithm would only require a 
number of steps polynomial in log(l/e). The origin of 
this limitation lies not only in the inability to design effi- 
cient measurements, but also in the accumulation of er- 
rors which occurs in the course of performing a quantum 
simulation. 

Here, we consider these limitations on the precision of 
results obtained by quantum simulations in the context of 
a specific algorithm for the simulation of pairing models, 
as proposed by Wu, Byrd and Lidar (WBL) [1] , which fol- 
lows the framework of the two methods described above. 
We present a study of the errors in its discrete time step 
implementation, and experimental results from a realiza- 
tion using a 3 qubit nuclear magnetic resonance (NMR) 
quantum computer, answering three questions: 1) What 
are the theoretical bounds on the precision of the quan- 
tum simulation? 2) How do faulty controls affect the 
accuracy of a simulation? 3) Can the theoretical bounds 
on precision be saturated by an NMR implementation? 

The WBL algorithm uses the classical FT algorithm 
described above to solve the question of the low-lying en- 
ergy gap in pairing Hamiltonians. Pairing Hamiltonians 
are used to describe both nuclear dynamics and super- 
conductivity [17-19] and are usually written in terms of 
Fermionic creation and annihilation operators and c 
as 

n N 

#pair = -j( c] rn C ™ +C -rn C -rn)+ ^ V ^ ( C ln C -m C l C -l) > 

m— 1 m,l—l 

where n is the total number of modes, e m is the onsite 
energy of a pair in mode m, and V m i are the coupling 
constants between modes. 

WBL map the pairing Hamiltonian onto the qubit 
Hamiltonian 

n v V 

#pair = ^ ~Y { — Z m ) + -^-(X m Xl + Y m Yi) , 

m=l m<l 

where X m ,Y m , and Z m are the Pauli operators on the 
rath qubit and v m = e m + Va (dropping an unimportant 
global energy shift, and using the standard convention 
Z\0) = |0)). The number of modes that can be sim- 
ulated equals the number of qubits n, and the number 
of pairs equals the total number of qubits in the state 
1). WBL show that for a specific number of pairs, one 
can approximately prepare the state \^i) by quasiadi- 
abatic evolution. Since i/ pa i r is a 2-body Hamiltonian, 
the system's evolution can be efficiently simulated on a 
quantum computer for any number of qubits [3]. WBL 
propose implementation of their algorithm using an NMR 
quantum computer, in which the operator M is simply Z 



for a single spin. An advantage of the ensemble nature of 
NMR is that a single measurement for a simulated time t 
yields (M(t)). Fixing a maximum energy width and de- 
sired precision makes the FT independent of the number 
of qubits. 

Let us begin by addressing the first question posed 
above, regarding theoretical bounds on the precision of 
this quantum simulation: how does the number of gates 
scale with the error e? The WBL method requires con- 
structing an operator that approximates the simulated 
Hamiltonian for times tk ■ The classical FT then yields an 
error of 2u E max /Q where E max is the largest detectable 
energy I /to (h = 1). In the case of using phase-estimation 
and the QFT, setting 2 q =Q yields the same precision. 
How long does it take to implement the Hamiltonian for 
a time Qto compared to implementing a Hamiltonian for 
time t$l In general, the operator is assumed to be con- 
structed of repetitions of the basic time step and requires 
Q more gates or time. This leads to the number of gates 
scaling inversely with the error. A similar problem faced 
in quantum factoring is overcome in Shor's algorithm by 
a clever way to perform the modular exponentiation [5] . 

A second bound on the number of gates required arises 
in calculating the time required to perform the algo- 
rithm. Quantum simulations typically employ a Trotter 
formula to approximate a Hamiltonian from combina- 
tions of non-commuting Hamiltonians [16]. For example, 
given the ability to evolve under Hamiltonians Ha and 
Hb, one can approximate evolution under Ha+Hb with 
bounded error. To lowest order, exp {—it{HA + Hb)) = 
(cxp(~itHA/k)exp(—itHB/k)) k + 5, where for 
\\[H Al HB}\\t 2 < 1, the error S is 0(t 2 /k). Higher- 
order techniques can yield an error 0{t m+1 / k m ) at the 
cost of needing 0(2 m ) more gates [20]. 

This approximation method leads to a subtle but im- 
portant difficulty in reducing the gate count for sim- 
ulations. It is apparent that the Trotter formula de- 
mands an exponential increase in the number of discrete 
gates for an exponential decrease in the error. However, 
from a Hamiltonian control perspective, this conclusion 
seems unfair, because the total time required can be 
small even if the gate count is high. Specifically, the 
gate UA(t/k) = cxp (—UHa/ k) requires 1/fc the time 
needed to implement J7a(*)- Therefore, the simple Trot- 
ter method given above requires only time 2t, indepen- 
dent of k. This implies that "Trotterization" errors in- 
volved in approximating desired Hamiltonians can be re- 
duced efficiently with respect to the time cost. 

Unfortunately, this optimistic observation is incompat- 
ible with fault tolerant error correction [21, 22], which 
will likely be needed to extend simulation times beyond 
limits imposed by qubit dccoherencc times. This is be- 
cause the fault-tolerant implementation of UA{t/k) takes 
approximately the same amount of time as the gate 
UA(t), whether using teleportation [23] or the Solvay- 
Kitaev approximation [24]. 
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Consequently, fault-tolerant simulations using the 
Trotter formula and the FT/QFT require a number of 
gates and amount of time that scales as 1/e 2 . Circum- 
venting this problem would require removing the ineffi- 
ciency of Trotterization, and constructing methods to ap- 
proximate Unit) with error e using poly(log(l/e)) gates. 
However, such methods would imply that the approxima- 
tion of Uuiqi) could take only poly (log (qj) more gates 
than the simulation of Unit)- Such a dramatic simplifi- 
cation may hold for specific H(t), but is unlikely to be 
possible for general H(t). The WBL algorithm studied 
here unfortunately does not scale efficiently when made 
fault tolerant. 

These theoretical bounds establish that present quan- 
tum simulations such as the WBL algorithm, using the 
QFT or the FT, require a number of gates which scales 
inversely with the desired answer precision for two rea- 
sons: fault-tolerant gate construction and the precision 
of a finite FT. Therefore, the time required for a <i-qubit 
quantum simulation is 0(n d /e r ), where r > 1 varies de- 
pending on the approximation methods employed, and 
r = 1 when quantum error correction and fault-tolerant 
gates are not used. 

We turn now to the second question, which concerns 
the impact of faulty controls in a real physical implemen- 
tation of the WBL algorithm. Recall that the foundation 
of the WBL algorithm is approximation of the unitary 
evolution under H pair , U pair (qt Q ) = cxp {-iH pa , iT qt )- An 
ideal NMR implementation accomplishes this by a re- 
peatable pulse sequence V pa ,i T (to), where U pa i r (qt ) ~ 
(V pa ,i r (to)) q . -ffpair contains three noncommuting parts: 
Ho = E m V -f(-Z m ), H XX = E^i^X^, and 
Hyy = J2 m <i ^'Ym.Yi. Assuming that the correspond- 
ing unitary operators Uo(t), Uxx(t), and Uyy(t) can be 
implemented, Vp a i r (io) can be constructed using the third 
order Trotter-Suzuki formula [20, 25] 

^pair(to) - [U (t /2k)U XX (t /2k)Uyy(t /k)x 

Uxx(h/2k)Uo(t /2k)} k , 
yielding an expected error \\U pa-ir (t ) - Vp air (io)|| = 

o(4/e). 

However, this ideal procedure is not actually achieved 
in a real experiment because the unitaries are not direct 
implementations of the Hamiltonians but instead com- 
posed from a series of pulses. These pulses depend on 
assumptions about the system Hamiltonian that become 
unreasonable for short simulated times. The reason is 
that all real systems have small, often unknown, energy 
shifts that are averaged away for large simulated times. 
When not using a fault-tolerant construction these shifts 
can lead to faulty controls. In atomic physics, for exam- 
ple, undesired Stark shifts need to be carefully accounted 
for in order to get exact rotations [26] . 

Control errors in NMR quantum computation arise, 
for example, since single qubit gates require finite time 
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FIG. 1: Three qubit quantum circuits for the unitaries Uxx 
(top) and Uyy (bottom), implemented using method Wl. 
These are depicted for simulating Hamiltonian Hi (see text). 
For Hi, -| pulses are applied to qubit c in parallel with those 
on a and b, and the decoupling X pulse is omitted. 



and unwanted two qubit coupling occurs during this 
time. In a static magnetic field B a z, the unitary evo- 
lution of a typical used spin system [11, 27] in the rotat- 
ing frame is given by Uzz(t) = cxp (-i J^ij ^JijZiZjtj , 
where the Jy are the scalar coupling constants. The 
time required for a radiofrcquency (RF) pulse to ro- 
tate individual spins by tt radians is much smaller than 
the typical delay times td during which no RF is ap- 
plied, ijr <C td ~ 1/Jij- Thus, it ordinarily suffices 
to approximate the RF pulses as (5-functions in time, 
implementing perfect single qubit rotations R^iO) = 
exp \^(Xi cos 4> + Yi sin </>)] . However, this approxima- 
tion breaks down as td becomes comparable to t^, causing 
the expected evolution to be best described not by dis- 
crete one and two-qubit gates, but instead by the piece- 
wise continuous time-dependent Hamiltonian H nmr (t) ~ 
9i{t)Xi + J2 t Mm + J2 iKj iJijZiZj. This discrep- 
ancy leads to additional errors in implementations of 
quantum algorithms and simulations, which, for a small 
number of qubits, can be mitigated using optimal control 
techniques[10, 28]. 

The impact of such control errors in an NMR im- 
plementation of the WBL algorithm can be stud- 
ied by comparing a baseline realization with no con- 
trol error compensation (denoted Wl) versus another 
with simple error compensation (denoted W2). The 
baseline Wl realization implements Uq using com- 
posite pulses to create rotations about the z axis, 
Uo = n m ^/ 2 (^/ 2 )- R o l (^» l )-R m ^/ 2 (^/2); an equivalent 
method, used elsewhere[l, 10] temporarily shifts the ro- 
tating frame. Control errors arise in the simulation of 
Uxx and Uyy, which are generated by applying single 
qubit pulses to rotate the scalar coupling from the z axis 
to the x and y axis, using the quantum circuit in Fig. 1. 
Control errors in this baseline realization are thus small 
only when delays needed to generate Uxx are long com- 
pared to the time required to perform single qubit gates, 
but also short enough that the Trotter error is small. 

A simple, scalable compensation technique for control 
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FIG. 2: Frequency-domain spectra of Hamiltonian H2 ob- 
tained using methods Wl (marked by circles) and W2 (di- 
amonds). Dots are experimental data, and solid lines are 
Fourier transforms of an exponentially decaying sinusoid 
(with four free parameters) fit to the time dependent data. 
The width of the exact curve is taken to be the dephasing 
rate (1/r) of the 13 C nucleus. 



errors provides a contrasting realization of the WBL algo- 
rithm for comparison. This W2 realization accounts for 
unwanted two qubit coupling during single qubit gates by 
reducing delay times during which coupling is desired. 
Specifically, every instance of R ( f >1 (9i)Uzz(t)R<f> 2 {0i) is 
replaced with R ( f >1 (9i)Uzz{t — a)R ( p 2 (9i), where a = 
+#2)- This technique was critical in the successful 
implementation of Shor's algorithm with NMR[29]; here, 
it is used with care, since many Hamiltonians have the 
same A as the pairing Hamiltonian of interest, and it 
is possible to tune a to get the right A for the wrong 
reasons. 

Numerical simulations comparing Wl and W2 show 
that the effect of such control errors on the WBL algo- 
rithm is a shift in the estimated gap value A from the ex- 
pected value. This shift can be quite significant, as shown 
in Fig. 2, and indeed can dominate errors due to other 
imperfections, such as the Fourier transform. Compen- 
sating for unwanted scalar couplings in NMR implemen- 
tations of quantum simulations is thus vital for obtaining 
correct results; implementations with other physical sys- 
tems will similarly have to deal with faulty controls. 

Finally, we consider the third question: saturation of 
the predicted precision bounds with an NMR implemen- 
tation of the WBL algorithm. The WBL algorithm is 
parameterized by the number of qubits, n, the simulated 
time step to, the number of steps Q, the degree of Trot- 
terization k, the adiabatic time step t a d, and the number 
of adiabatic steps S. An NMR system is characterized 
by a characteristic decoherence time r, and following the 
discussion above, it is convenient to work with a small 
number of qubits for times shorter than ~ 3r, such that 
quantum error correction is unnecessary. We implement 



two specific instances of the pairing Hamiltonian £f pa i r 
involving three modes (n = 3) and two pairs, leading to 
a 3-dimensional Hilbert space spanned by 1 101} , 1 110) , 
and |011). The simulation is started the ground state of 
the two spin up subspace of Ho, |011), prepared using 
temporal labeling [30]. 

WBL estimated the expected size of the system that 
could be simulated without error correction by choosing 
k/t = 0.1A and eft — A. They found that the number 
of gates required scales as in^A/epT, including the nec- 
essary decoupling pulses. The gate time, t g is assumed 
equal to 10~ 5 r and for up to n = 10, A can be found to 
precision e w A. Here, we find A to precision e « A/100, 
and the number of qubits is n = 3, consistent with the 
WBL bound n < 4 for these parameter choices. 

The first stage of the WBL algorithm is to quasiadia- 
batically evolve into the ground state of -ff pa i r , with 
discrete changes in the simulated Hamiltonian, using a 
procedure previously demonstrated [31]. The Hamil- 
tonian used at each discrete step s is H a d(s) = (1 — 
s/S)H + {s/S)H pa i I , where S is the maximum number 
of steps. Unitary evolution U ac i at each step for time t a d 
is then approximated using the above pulse sequences, 
as U ad = ]lf=o v ad(s, tad)- Preparation of the state 
requires evolving at a rate faster than that for adiabatic 
evolution, thereby exciting the state \Ei). This quasia- 
diabatic evolution is accomplished by reducing £ or t a d 
compared to the adiabatic case [1] . Higher-energy states 
will also be excited, but S and t a d can be adjusted to 
minimize this. Quasiadiabatic evolution in this experi- 
ment was attained with S ~ 4 steps and t a d = 1/700 
s. Note that for \\H X x + H Y y\\ > \\H \\ there can be 
a phase transition as s is changed [32]; as the gap goes 
to zero at the phase transition this can be problematic, 
since the number of steps required for successful quasia- 
diabatic evolution grows inversely with the gap. 

The second stage of the algorithm is evolution of the 
state \^i) under the pairing Hamiltonian for Q timesteps 
of duration t with k = 2. These parameters are chosen 
such that l/{Qt ) w A/100 , Qt < t, and k/t > 0.1A. 
Note that many Q and to yield the same tpr', this is used 
to our advantage below. 

We performed our experiments using a 500 MHz Var- 
ian vmTY INOVA spectrometer and 13 C-labeled CHFBr 2 , 
with coupling strengths Jhc = 224 Hz, Jhf = 50 
Hz, and Jqf = —311 Hz. The two pairing Hamiltoni- 
ans simulated were Hi, the "natural" Hamiltonian, in 
which V12 = ttJhc, V13 = kJ hf , and V23 = ttJcf, 
and a harder case, an artificially constructed Hamilto- 
nian Hi, in which V a b = ttJhc and V ac = Vbc = 0. For 
both Hamiltonians, vi = 1507T Hz, u 2 = 1007T Hz, and 
z/3 = 507T Hz. Implementation of Hi required an addi- 
tional decoupling pulse compared to Hi, following the 
circuit in Fig. 1. 

Each Hamiltonian was simulated for times to to QU), 
using the Wl and W2 methods, and an NMR spectrum 
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Model 


A/Hz 


Method 


A exp /2n-}lz 


T e / HIS 


io/ms 


Q 


Hi 


218 -27T 


Wl 


227 ± 2 


180 


1 


400 






Wl 


220 ±2 


250 


2 


200 


H 2 


452 -27T 


Wl 


554 ± 10 


30 


.5 


200 






W2 


440 ±5 


80 


.5 


200 



TABLE I: Experimental results for gaps found for Hamiltoni- 
ans Hi and Hi. Estimated gaps (A exp ) and effective coher- 
ence times (r e ) for given time steps to and number of steps 
Q are obtained by least-squares fitting of the time-dependent 
NMR peaks to an exponentially-decaying sinusoid. 



was acquired for each time duration. A classical discrete 
FT of the NMR peak intensities of one spin (hydrogen) 
over the t& yielded four spectra of i? pa ir- The experi- 
mental result A exp was determined by a least-squares fit 
of the highest-signal NMR peak to a damped sinusoidal 
function with frequency A and decay rate r e . 

Ideally, the result should find A = 218 • 2?r Hz for H u 
and A = 452 • 2ir Hz for H 2l as determined by direct 
diagonalization. Note that for H 2 , A is the energy differ- 
ence between \G) and I-E2}, since \Ei) is not connected 
by usual adiabatic evolution; the larger gap requires that 
to/kbe smaller when simulating H 2 ■ For the experimen- 
tal result, we expect that A exp —A + e sys ± cft where 
e sys is an offset due to Trottcrization and/or faulty con- 
trols. Q and to determine the theoretical bound on the 
precision, in the absence of control errors, and the exper- 
iment should saturate this bound when the A exp is cft 
from the actual value A. 

Experimental results for the spectra of H 2 are shown 
in Fig. 2; experimental parameters Q and to and numer- 
ical results from the analysis for A exp and r e for each 
experiment are summarized in Table I. 

The impact of systematic and random errors was in- 
vestigated by simulating Hi with Wl (no control error 
compensation) for cft = 2.5 • 27r Hz at two different sim- 
ulation times, t = 1 ms an d to = 2 ms. As expected, 
the random error for both cases is « cft- Note that the 
systematic error increases with smaller to- This signals 
that the error due to unwanted scalar coupling becomes 
larger than the errors due to the Trotter approximation. 
Consequently, a slightly longer to yields a systematic er- 
ror that is within cft of the exact answer, saturating the 
predicted theoretical bounds on precision. 

Convergence to the correct result is another important 
issue for all discrete time simulations. For this 3-qubit 
system, we performed a detailed numerical simulation to 
determine that to = 2 was optimal. For a large system 
this is no longer possible, and convergence tests would 
need to be used to verify the answer. The procedure 
would reduce to ( OT increase k) until the change in A exp 
was smaller than the desired precision. 

While the results for Hamiltonian Hi were good even 



without control error compensation, the effects of control 
errors were very evident in the results for H 2 . Hamilto- 
nian H 2 was implemented with Wl (no error compen- 
sation) and W2 (simple error compensation) for epT — 
10 • 2-7T Hz and t = 0.5 ms. The shorter time step was 
necessary because the larger A made the simulation more 
sensitive to Trotter errors. Comparing the W2 and Wl 
results shows that with no control error compensation, a 
gap A is found that is A/5 away from the actual value. 
In contrast, with simple error compensation A is epT 
from the actual value, saturating the theoretical bound. 
Future implementations should certainly strive to detect 
and bound control errors; this could be done by verifying 
that A exp scales as tp for small values of to, as theoreti- 
cally expected. 

In conclusion, we have studied the theoretical and em- 
pirical bounds on the precision of results obtained with 
quantum simulations, in the context of the pairing Hamil- 
tonian algorithm proposed by Wu, Byrd, and Lidar. We 
have implemented the smallest problem instance that re- 
quires quasiadiabatic evolution, verifying that the algo- 
rithm computes the gap A to within the precision of the 
method. We also find, however, that simulations of this 
type are particularly sensitive to systematic errors in the 
applied Hamiltonian and that fault-tolerant implementa- 
tions are inefficient with respect to precision using cur- 
rent Trotter approximation methods. 

Nevertheless, in practice, when only limited precision is 
desired and for a sufficiently large system, quantum sim- 
ulations may still outperform classical numerical simula- 
tion, as demonstrated for molecular energies [33] . Avoid- 
ing the cost of precision is desirable, and can be done by 
designing quantum simulations to explore questions that 
are insensitive to the microscopic details of the Hamil- 
tonian [34]. How to develop quantum simulations for 
faulty small scale (10-20 qubit) quantum computers that 
can outperform classical computations remains an open 
question. 
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